Correlation Dimension of Inertial Particles in Random Flows 
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We obtain an implicit equation for the correlation dimension D2 of dynamical systems in terms 
of an integral over a propagator. We illustrate the utility of this approach by evaluating D2 for 
inertial particles suspended in a random flow. In the limit where the correlation time of the flow 
field approaches zero, taking the short-time limit of the propagator enables D2 to be determined 
from the solution of a partial differential equation. We develop the solution as a power series in a 
dimensionless parameter which represents the strength of inertial effects. 
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The behaviour of small particles moving independently 
in complex flows is a fundamental problem in fluid me- 
chanics, which has applications in understanding rainfall 
[H , planet formation [2, Q and many areas of technology 
and environmental science. It is known that when the in- 
ertia of the particles is significant, clustering may occur 
Q , which can lead to an increase in the rate of collision or 
aggregation of the particles, and which can also affect the 
scattering of electromagnetic radiation. In developing a 
description of these processes the most natural way to 
quantify the clustering is to consider the number of par- 
ticles M inside a ball of radius 5r centred on any given 
particle. If this quantity has a power-law dependence for 
small 6r of the form M ~ Sr^'^ (with D2 less than the 
dimension of space, d), the particles cluster onto a frac- 
tal attractor. The quantity D2 is termed the correlation 
dimension Q. The clusteringprocess is in fact found to 
approach a fractal attractor [6| . 

It is desirable to develop a theoretical understanding 
of the clustering effect. It has been ascribed to parti- 
cles (which we assume to be much denser than the fluid) 
being centrifuged away from vortices Q, but other ex- 
planations (for example, caustics 0,01) are possible. In 
particular, a model with a short-time correlated velocity 
field, analysed in P, gives good agreement with a nu- 
merical determination of the Lyapunov dimension of 
particles in Navier-Stokes turbulent flow, reported in 
(the Lyapunov dimension was introduced in and is 
discussed in fsj). The task of calculating the more physi- 
cally interesting dimension D2 by analytical methods has 
appeared to be intractable, but we show that D2 can be 
obtained more easily than Dl- We give a general pre- 
scription for calculating the correlation dimension, which 
can also be applied to other types of dynamical system. 
We show that when the turbulent velocity is modelled 
by a random vector field with a short correlation time 
(that is, for the model analysed in Q), this leads to an 
expansion of D2 as a power series in a parameter e which 
is a dimensionless measure of the inertia of the particles. 
The coefficients of this series may be obtained exactly to 
arbitrarily high order. We show how convergent results 
are obtained using a conformal Borel summation. 

The correlation dimension D2 may be defined in terms 



of the expected number M{5r) of particles inside a ball 
of radius 5r surrounding a test particle: 



D2 = lim 



ln[(AA(fr))] 
\n{6r) 



(1) 



(where {X) denotes an average of X), provided this limit 
exists and satisfies D2 < d, where d is the dimensionality 
of space. This implies that {M{5r)) ~ Sr^'^ which is the 
radial part of the volume element of a ball in D2 dimen- 
sions. If the limit in ([T]) is greater than or equal to d, there 
is no clustering, and D2 = d. While D2 has fundamen- 
tal importance, it is difficult to calculate analytically. It 
can be expressed in terms of the large deviation statistics 
of the finite-time Lyapunov exponents, a{t) [13, [T^ - 
These statistics are very difficult to calculate by means 
other than numerical simulations (although they have 
been evaluated for the Kraichnan model for advection 
in short-time correlated flows [IBl)- Earlier studies of D2 
for particles with significant inertia have been numerical 
evaluations [l3.[l5|. 

We consider the motion of small, dense particles sus- 
pended in a turbulent fluid with velocity field u{r,t). 
The motion of a particle at position r moving with ve- 
locity V is determined by viscous damping of the particle 
relative to the fluid. The equations of motion are 



V = —^[v — u{r{t), t)] 



(2) 



where we use the notation X = dX/dt and where 7 is a 
damping rate proportional to the viscosity. In this paper 
we consider how to extract information about D2 from 
a quantity Zi(t) which is defined to be the logarithmic 
derivative of the separation Sr between two particles: 



or 



(3) 



An equation of motion for Zi which is valid when 6r is 
sufficiently small may be obtained from the linearisation 
of ^ as discussed below: Zi{t) may be coupled to one 
or more additional variables Z2{t), . . ., but the equations 
for the Zi are independent of Sr provided that quantity 
is sufficiently small. We also consider the variable 



Yit) = Infr(t) 



(4) 
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which is related to Zi hy Y = Zi. Note that Y is re- 
lated to the finite-time Lyapunov exponent a{t) at time 
t: we have Y{t) — Y{0) = ta{t) (provided 5r is everywhere 
sufficiently small). In the following we discuss the two- 
dimensional case where Zi is coupled to one additional 
variable Z2. We consider the joint probability density 
p{Y, Zi, Z2) of y, Zi and Z2. Because the equation of 
motion of Zi and Z2 is independent of y = In 5r when 
the linearised equation is valid, in the steady state the 
joint distribution factorises, with the distribution of Y 
being in a form which reflects the translational invari- 
ance in Y . Because the eigenfunctions of translations are 
exponential functions, the steady-state joint distribution 
of r, Zi, Z2 is 



p(r, Zi, Z2) - exp(ar)pz(Zi, Z2) 



(5) 



for some constant a. This form is not normalisable, but 
it should be remembered that ([5]) is only valid when 
Sr is sufficiently small. In the case where a > 0, 
the form ([5]) can be matched to a distribution which 
is valid for large 5r to make a normalisable solution, 
whereas a < is not allowed. The distribution ([5]) im- 
plies that the distribution of Y has probability element 
dP = exp(Q;F)dF — 5r°'~^d5r. The relation ([T|) implies 
that the probability for the separation to be in an interval 



A5r is dP = 5\ 



,-D2-l 



d(5r, so that 



D') = a 



(6) 



The condition for determining D2 = a is that this dis- 
tribution ([5]) should be invariant under time evolution. 
This is expressed in terms of a propagator for the time- 
evolution of Y and Z = (Zi, Z2). Specifically, this prop- 
agator if (AY, Z, Z' , At) is defined to be the probability 
density for Y to change by AF and for Z = (^1,^2) 
to change from Z' to Z in time At. Stationarity of the 
distribution ^ then leads to 



Pz{Zi,Z2 



dAy 



dZ[ 



dZ'2 



X exp(-«Ay) K{^Y, Z, Z\ At) pz{Z[,Z^) (7) 

which is satisfied for all At. In the case At 00, the 
propagator K is related to the large-deviation probability 
density function for the finite-time Lyapunov exponent. 
This leads to a formulation (to be discussed in a later 
paper) which is equivalent to some earlier theories for de- 
termining D2 fsl. [l^. [ist . Here, however, we concentrate 
upon the short-time limit. At 0. We shall see that 
this leads to an analysis of D2 in terms of a differential 
equation, which is much more analytically tractable. 

To make further progress we need to consider the 
equation of motion for the variables Zi , Z2 in the two- 
dimensional case. Parts of the calculation follow [l^ , but 
here we use a simpler operator algebra. The linearised 
equations of motion corresponding to are Sr = Sv and 
Si) = —jSv + 7E(t)i5r where E(t) is a 2 x 2 matrix with 



elements Eij{t) = dui/ drj{r{t),t). We write Sr = Srng 
and Sv — ZiSrxig + Z2(5rng_|_^/2j where is unit vec- 
tor in direction 0. Expressing the linearised equations of 
motion in terms of the variables Sr, Z\, Z2 we obtain (l6l 



Z\ = -7Z1 + (Z2 - J 
Z2 = ^7^2 — 2Z1Z2 



(8) 



where Ed{t) = ng • E(t)n0 and Eo{t) = ne+^/2 • E(t)ne, 
and Sf = ZiSr, 6 = Z2 (so that the definition of Zi is 
consistent with ([3])). It might be expected that the dis- 
tribution of (Zi, Z2) obtained from the long-time limit of 
the evolution of equation ([8]), which we term pq{Zi, Z2), 
is the same as the distribution pz{Zi, Z2) in ([7]). How- 
ever, pz differs from po because it is conditioned upon 
being at a particular value of Y. If a > 0, particles 
reaching a negative value of Zi arrive from a larger value 
of y, where the probability, density is larger. This im- 
plies that the distributions po and pz are different, and 
that Pz has a smaller mean value of Zi than po- 

Next we must specify a model for the two-dimensional 
velocity field u{r,t). We allow this to be partially com- 
pressible by writing u = V$ + V A ^'63. In order to 
use statistical techniques we consider the stream func- 
tion ^'(r,t) and potential $(r,t) to be random scalar 
fields with specified correlation functions. We shall as- 
sume that {<i>{r,t)<i>{r',t')) = C{\r - r'\, \t - t'|), where 
C{R,t) has support ^ (the correlation length) and t (the 
correlation time) in R and t respectively. Also, we as- 
sume that and 5* are uncorrelated and that the cor- 
relation function of \1/ is proportional to that of such 
that (\E'^)/($^) = for some number f3. Furthermore, 
in this paper we consider the limit where we the correla- 
tion time T is sufficiently small that the randomly fluc- 
tuating terms in ([5]), Ed{t) and Eo{t), can be treated as 
white noise. In this case the equations of motion for Zi, 
Z2 become a pair of coupled Langevin equations, and the 
probability density po(^i, ^2) generated by equation ^ 
obeys a diffusion equation, which can be written formally 



dpo ^ 

where is a Fokker-Planck operator: 



(9) 



•^opo - ^[(7^i + ^?-^|)Po]+Pii^ 

+ ^^[(^Z^ + '^Z^Z^)P^\+T^^^^ ■ (10) 

Here the diffusion coefficients are expressed in terms of 
correlation functions of the velocity gradients: 



/oo 
dt {Ea{t)E,im 
-00 



(11) 



Now we consider how equations ([H]), (jlOp are 
used to construct the short-time propagator in ([7]). 
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For small At, Y evolves ballistically, with veloc- 
ity Zi Z[. In the short time limit, the 
action of the propagator K {AY, Z , Z' , At) in ([8|) 
on a function f{Y,Zi,Z2) can therefore be writ- 
ten as fK{Y,Z^,Z2) = f{Y - Z^At,Z^,Z2) + 
AtTaf{Y,Zi,Z2) + 0{At^). The equation (0) de- 
termining self-reproduction of pz{Zi,Z2) therefore 
becomes pz{Zi,Z2) — exp{-aZiAt)pz{Zi, Z2) + 
AtTopziZi,Z2) +0(Ai2). Extracting the 0{At) term 
gives the differential equation 



aZipz{Zi, Z2) - Tqpz(Z\, Z2) = 



(12) 



Upon integrating over space, and using the fact that the 
operator J-q is a divergence, we have 



/oo />oo 
dZi / dZ2 Z,pz{Z^,Z2) = (Z,) 
-00 J — 00 







(13) 



The value of D2 is determined by finding the value of 
a for which a normalisable solution of (|12|) can be ob- 
tained for which the mean value of Zi is zero. The equa- 
tions (fT2)) and together constitute an exact method 
for determining D2 — a. Their extension to the three- 
dimensional case is straightforward. 

It is useful to make a change of variable from {Zi, Z2) 
to scaled variables {xi,X2) defined by Xi = \J ^ jVaZi, 
and to use a dimensionless time t' = 'jt. We also intro- 
duce two dimensionless parameters, e, which measures 
the importance of inertial effects, and F, which is a con- 
venient measure of the relative magnitudes of 4" and <&: 




F = 



22 



1 + 3(3^ 
3 + /32 



(14) 



Using these new variables ([T^ becomes an equation for 
the joint probability density P{xi,X2) of xi, X2'. 

FP^Q^^^[{xi+e{xl-Txl))P] 



d 

+ i^[(a;2 +2e.xiX2)P] + 
0x2 



Q2p Q2p 



dx{ 



dxl 



eaxiP (15) 



(which defines the differential operator F{e, a, F)). Equa- 
tion (|15p is to be solved with the supplementary condition 
(xi) = 0, which can only be satisfied for isolated values 
of a. Our solution below obtains one unique value of a, 
which is D2. 

We now develop the solution as a series expansion in 
e, using a system of annihilation and creation operators 
which are analogous to those used in quantum mechan- 
ics. We use a notation similar to the Dirac notation, 
whereby a function f{xi,X2) is denoted by a vector |/). 
We expand both the solution \P) of and the value 
of a for which the solution of this equation exists and 
satisfies (xi) = as power series in e: 



fc=0 



E 

fc=0 



(16) 



We write the Fokker-Planck operator in ([151) ^ 
F ^ Fo + €{G - axi) 



(17) 



(thereby defining operators Fq, G). The unperturbed 
steady-state |Po) satisfying -Fbl^o) = is Pq{xi,X2) — 
exp[— (a;^ +x\)/2]/2n, and other eigenfunctions of Fq are 
generated by creation operators di and annihilation op- 
erators 5,-: 



(18) 



These operators generate eigenfunctions satisfying 
FQ\4>nm) — ~{n + m)\4>nm.), accordiug to the rules 

ai\4'n,rn) = \4>n+l,rn) ^l|'/'n,m) = '^|0n^l,m) 
a2\4'n,rn) = \4'n,rn+l) b2\(l)n,m) = 'm\(l)n,m-l) (19) 

with l^oo) — \Po)^ which is normalised as a probability 
density. The states \Pk) in ([16]) with be expressed as 
linear combinations of the eigenfunctions | (/)„,«): 



oo oo 



(20) 



n— m— 



In general these eigenfunctions are neither normalised, 
nor do they form an orthogonal set, but these properties 
are not required in the following arguments. We first con- 
sider the implications of the requirement that {xi) — 0. 
Using (|18p and (|19p . by an inductive argument involving 
repeated integration by parts we have: 



dxi 



<^X2 4>nm{xi,X2)Xi = (S„l(5mO (21) 



so that the condition {xi) = is satisfied by requiring 

that pf^ = in ([20|) for all k. 

Substituting ([TO)) into ([T5)) leads to a recursion giving 
\Pn) in terms of all of the preceding approximations: the 
term of order e" is 

= Fo|P«) + [G - ao(ai + Si)] |P„_i) . . . 

- aj(ai+6i)|P„_i_j)...-a„_i(ai + 6i)|Po) -(22) 

There are two unknowns in this equation, |P„) and a„_i; 
all of the other |P, ) and aj are assumed to have been de- 
termined at previous iterations. For any value of an-i, 
equations ([22]) can be solved formally for |P„) by mul- 
tiplying by Fq^. For a state \Q) with coefficients qnm. 

we have Pq^^IQ) = - ^^=0 Em=0 TITT^T^nml'/'nm)- The 

action of Fq^ upon a general state \Q) is therefore unde- 
fined unless the coefficient goo is equal to zero. At each 
order we can solve ([^^ for |P„) choosing the value of 

a„_i so that p^g^ = 0. Note that the operator G con- 
tains raising operators as left factors, so that Pq"^G|/) 
exists for any state |/). However, because there is a low- 
ering operator bi acting on the states |Pfc), the action of 
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FIG. 1: Correlation dimension D2 of the model ([2]), as a func- 
tion of the dimensionless measure of inertia e, defined by (I13|l . 
Here F = 3 (incompressible flow) and r ^ (rapidly fluctu- 
ating flow field). The numerical data (o) are compared to 
the quadratic approximation of (|23|) . (dotted line) and to the 
Borel summation of the series (|22|l . (solid line, red online). 
The summation was evaluated using the method in , usmg 
the conformal map z — 2u''/s(l — m)", with 1/ — f/4, s — 25. 
The results converge as the number of terms used in the Borel 
summation, fcmax, increases: the curves for fcniax = 10, 20, 30, 
40, 50 lie on top of each other. 

multiplying the terms in by F^^ is only defined if all 

of the \Pk) are chosen so that p[''q = 0. However, we have 
already seen that this is precisely the condition to ensure 
that the solution satisfies (xi) — 0, that is, the solvabil- 
ity condition upon (22) coincides with the condition 
The generation of the series was automated using an 
algebraic manipulation program. Iterating the equation 
using the initial condition \Po) = \(j)oo) leads to the 
following series expansion for 1)2 (e): 

(23) 

All aj with odd j are equal to zero, and all the co- 



efficients are zero when F = 1. For F = 3 (so that 
V • M = 0) the first few non- vanishing coefficients are 2, 
-24, 528, -28800, 1654848, -128860416, so that the se- 
ries is clearly divergent with alternating signs. It is inter- 
esting to consider whether this series contains a complete 
description of D2{e)- We investigated its evaluatio n by 
means of a Borel summation technique described in [17| . 
The Borel transform B{z) = J2'kLo('^k/kl)z'' of D2{e) is 
convergent inside a disc (of radius 1/12), but inversion 
of B{z) to yield D2{e) requires its Laplace transform, 
which is an integral over z G {0,oo). This is facilitated 
by making a conformal transformation to a new vari- 
able u, defined by z = 2'^u/s{l — u)" (where s are 
constants), so that the positive z axis is mapped to the 
interval u S (0,1). We find that the expansion of B{z) 
as a series in u has decreasing coefficients when v = j 
and s = 25 (indicating that B{z) is analytic in the image 
of the disc |u| < 1). Performing the integral in the u 
variable gives a summation of the series which converged 
as the number of terms, /cmax, was increased. Figure [1] 
illustrates the results for F = 3. For small e there is ex- 
cellent convergence to a numerical evaluation of D2{e). 
For large e, however, while the Borel summation con- 
verges as fcmax is increased, it diverges from the numeri- 
cal evaluation. This indicates that there is a component 
of D2 (e) which has no representation as an analytic func- 
tion. Non-perturbative approaches to equation (|12p are 
required to describe this non-analytic contribution. 

Equation ([7]) can be used to determine the correlation 
dimension of other stochastic dynamical systems, includ- 
ing cases where the random component has a finite cor- 
relation time, and also for deterministic systems. A full 
account of the use of equation ([7]) to determine the cor- 
relation dimension will be published elsewhere. 
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